rm(list = ls())
library(DropletUtils)
## Warning: package 'DropletUtils' was built under R version 4.1.2
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## Warning: package 'matrixStats' was built under R version 4.1.2
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: GenomicRanges
## Warning: package 'GenomicRanges' was built under R version 4.1.2
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Warning: package 'GenomeInfoDb' was built under R version 4.1.2
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
library(scater)
## Loading required package: scuttle
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.1.2
library(scran)
## Warning: package 'scran' was built under R version 4.1.2
library(ggplot2)
library(Seurat)
## Warning: package 'Seurat' was built under R version 4.1.2
## Attaching SeuratObject
## 
## Attaching package: 'Seurat'
## The following object is masked from 'package:SummarizedExperiment':
## 
##     Assays
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.2
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:Biobase':
## 
##     combine
## The following objects are masked from 'package:GenomicRanges':
## 
##     intersect, setdiff, union
## The following object is masked from 'package:GenomeInfoDb':
## 
##     intersect
## The following objects are masked from 'package:IRanges':
## 
##     collapse, desc, intersect, setdiff, slice, union
## The following objects are masked from 'package:S4Vectors':
## 
##     first, intersect, rename, setdiff, setequal, union
## The following objects are masked from 'package:BiocGenerics':
## 
##     combine, intersect, setdiff, union
## The following object is masked from 'package:matrixStats':
## 
##     count
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(viridis)
## Warning: package 'viridis' was built under R version 4.1.2
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.1.2
library(EnhancedVolcano)
## Loading required package: ggrepel
## Warning: package 'ggrepel' was built under R version 4.1.2
## Registered S3 methods overwritten by 'ggalt':
##   method                  from   
##   grid.draw.absoluteGrob  ggplot2
##   grobHeight.absoluteGrob ggplot2
##   grobWidth.absoluteGrob  ggplot2
##   grobX.absoluteGrob      ggplot2
##   grobY.absoluteGrob      ggplot2
library(formattable)
## 
## Attaching package: 'formattable'
## The following object is masked from 'package:scater':
## 
##     normalize
## The following object is masked from 'package:BiocGenerics':
## 
##     normalize
library(pheatmap)
library(ape)
## Warning: package 'ape' was built under R version 4.1.2
## 
## Attaching package: 'ape'
## The following object is masked from 'package:dplyr':
## 
##     where
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##   genomic data. Bioinformatics 2016.
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
##    in the original pheatmap() are identically supported in the new function. You 
##    can still use the original function by explicitly calling pheatmap::pheatmap().
## 
## Attaching package: 'ComplexHeatmap'
## The following object is masked from 'package:pheatmap':
## 
##     pheatmap
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.1.2
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:formattable':
## 
##     area
library(scales)
## Warning: package 'scales' was built under R version 4.1.2
## 
## Attaching package: 'scales'
## The following objects are masked from 'package:formattable':
## 
##     comma, percent, scientific
## The following object is masked from 'package:viridis':
## 
##     viridis_pal
library(wesanderson)
library(ggrepel)
library(ggeasy)
## Warning: package 'ggeasy' was built under R version 4.1.2
library(data.table)
## Warning: package 'data.table' was built under R version 4.1.2
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following object is masked from 'package:SummarizedExperiment':
## 
##     shift
## The following object is masked from 'package:GenomicRanges':
## 
##     shift
## The following object is masked from 'package:IRanges':
## 
##     shift
## The following objects are masked from 'package:S4Vectors':
## 
##     first, second
library(scuttle)
library(SingleR)
## Warning: package 'SingleR' was built under R version 4.1.2
library("Seurat")
library("monocle3")
## 
## Attaching package: 'monocle3'
## The following objects are masked from 'package:Biobase':
## 
##     exprs, fData, fData<-, pData, pData<-
library("org.Hs.eg.db")
## Loading required package: AnnotationDbi
## Warning: package 'AnnotationDbi' was built under R version 4.1.2
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
library("dplyr")
addTaskCallback(function(...) {set.seed(100);TRUE})
## 2 
## 2

#Load object

combined_sce1 <- readRDS("FinalObj.RDS")
combined_sce2 <- combined_sce1
sce2 <- as.Seurat(combined_sce2)
Idents(sce2) <- sce2$clusters
#
sce1 <- sce2
sce1
## An object of class Seurat 
## 31719 features across 43436 samples within 1 assay 
## Active assay: originalexp (31719 features, 0 variable features)
##  2 dimensional reductions calculated: PCA, UMAP

#UMAPs For Figure 5A and 5B

color_samp <- wes_palette("Darjeeling1",3)
wt_col <- color_samp[2]
ko_col <- color_samp[3]
my_color_palette <- hue_pal()(length(levels(sce1$clusters)))
DimPlot(sce1) + ggtitle("UMAP by Clusters")

DimPlot(sce1, group.by = "condition", cols = c(wt_col, ko_col), label = FALSE) + ggtitle("UMAP by WT vs KO") + theme(legend.text=element_text(size=13))

#MAKE HEAT MAP, CLUSTERPLOTS, VLN

markers <- FindAllMarkers(sce1, min.pct = 0.15, min.diff.pct = 0.04, logfc.threshold = .25, test.use = "MAST") #Try with #MAST
## Calculating cluster 1
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster 2
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster 3
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster 4
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster 5
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster 6
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster 7
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
top20 <- markers %>%  group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
sce1 <- ScaleData(object = sce1, block.size=100)
## Centering and scaling data matrix
temp3 <- GetAssayData(sce1, slot = "scale.data")
mat<- temp3[unique(c("JUN", top20$gene)), ] %>% as.matrix()
heatmapanno <- HeatmapAnnotation(Group = sce1$condition, Cluster = Idents(sce1), col = list(Group=c("WT" = wt_col, "KO" = ko_col), Cluster =c("1" = "#F8766D", "2" =  "#C49A00", "3" =  "#53B400", "4" =  "#00C094", "5" =  "#00B6EB", "6" =  "#A58AFF", "7" ="#FB61D7")))
col_fun = circlize::colorRamp2(c(-1, 0, 2), c("#140b34" , "#84206b", "#f6d746"))

##Figure 5D

Heatmap(mat, name = "Scaled Expression",  
        column_split = Idents(sce1),
        #column_title = "Cluster",
        cluster_columns = FALSE, 
        show_column_dend = FALSE,
        cluster_column_slices = FALSE,
        column_title_gp = gpar(fontsize = 7),
        column_gap = unit(.6, "mm"),
        cluster_rows = TRUE,
        show_row_dend = TRUE,
        row_dend_width = unit(20, "mm"),
        col = col_fun,
        row_names_gp = gpar(fontsize = 8),
        column_title_rot = 0,
        top_annotation = heatmapanno,
        show_column_names = FALSE,
        use_raster = TRUE,
        raster_quality = 4)

#Figure 5F sce3 <- subset(sce1, ident = 7, invert = TRUE) v1 <- VlnPlot(sce3, “BAX”, split.by = “condition”, split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0) v2 <- VlnPlot(sce3, “TNFRSF12A”, split.by = “condition”, split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0) v3 <- VlnPlot(sce3, “JUN”, split.by = “condition”, split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0) v1+v2 + v3

##Figure 5E

Idents(sce1) <-sce1$condition 
markers1 <- FindAllMarkers(sce1, min.pct = 0.05, min.diff.pct = 0.02, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster KO
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
top <- markers1 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
  top$avg_log2FC < -.10, '#F2AD00',
  ifelse(top$avg_log2FC > .10, '#00A08A',
         'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
apotosis1 <-c("FAS",
              
              "TNFRSF10A",
              "TNFRSF10B",
              "TNFRSF10C",
              "TNFRSF10D",
              "TNFRSF11B",
              "TNFSF10",
              "TNFRSF1A",
              "TNFRSF12A",
              "FADD",
              "CFLAR",
              "CASP1",
              "CASP2",
              "CASP3",
              "CASP4",
              
              "CASP6",
              "CASP7",
              "CASP8",
              "CASP9",
              "CASP10",
              "CASP14",
              "NAIP",
              "BIRC2",
              "BIRC3",
              "XIAP",
              "BIRC5",
              "BIRC6",
              "BIRC7",
              "BCL2",
              "MCL1",
              "BCL2L1",
              "BCL2L2",
              "BCL2A1",
              "BCL2L10",
              "BAX",
              "BAK1",
              "BOK",
              "BID",
              "BCL2L11",
              "BMF",
              "BAD",
              "BIK",
              "HRK",
              "PMAIP1",
              "BNIP3",
              "BNIP3L",
              "BCL2L14",
              "BBC3",
              "BCL2L12",
              "BCL2L13",
              "APAF1",
              "CYCS",
              "DIABLO",
              "HTRA2",
              "AIFM1",
              "ENDOG",
              "CARD8",
              "CARD6",
              "NOX5",
              "TP53",
              "CDKN1A",
              "CDKN1B",
              "CDKN2A",
              "CDKN2B")
ofinterest <- c("MKI67",
                "TOP2A",
                "NES",
                "MAP2",
                "TH",
                "EN1",
                "NR4A2",
                "HIF1A",
                "TP53",
                "RELA",
                "NFKB1",
                "NFKB2",
                "JUN",
                "FOS", "LMO3", "HES5")
geneofinterest <- c(ofinterest, apotosis1)
samplabels1 <- c("PEG10", "DYNC1I2", "JUN",
                 
                 "NMT1",
                 "PPDPF",
                 "PHPT1",
                 "NDUFA13",
                 "NDUFC2",
                 "NDUFA11",
                 "NME2",
                 "NDUFA12",
                 "RPS27L",
                 "INSM1",
                 "FABP5",
                 "HSPA1A",
                 "WLS",
                 "MICOS13",
                 "ASCL1",
                 "SOX2")
sce4 <- subset(sce1, ident == 3 | ident == 5 | ident == 6)
Idents(sce4) <-sce4$condition
markers2 <- FindAllMarkers(sce4, min.pct = 0.01, min.diff.pct = 0.01, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster KO
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
top2 <- markers2 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
  top2$avg_log2FC > .10 & top2$p_val_adj < .01, '#00A08A',
  ifelse(top2$avg_log2FC < -.10 & top2$p_val_adj < .01, '#F2AD00',
         'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
samplabels<- c("RPS27L", "PHLDA3", "TUBB2A","TUBB2B","RPL27A","PHPT1","CDKN1A","BAX",
      "NMT1",
      "BBC3",
      "CDKN2B",
      "TNFRSF10B")
EnhancedVolcano(top2,
                lab = top2$gene,
                x = 'avg_log2FC',
                y = 'p_val_adj',
                title = 'WT vs KO Top DE Genes',
                pCutoff = .01,
                FCcutoff = .10,
                pointSize = 1.0,
                labSize = 3.3,
                selectLab = samplabels,
                col=c('blue', 'green', 'black', 'red3', "yellow", "brown", "grey"),
                colAlpha = .8,
                colCustom = keyvals,
                drawConnectors = TRUE,
                widthConnectors = 0.1,
                lengthConnectors = unit(0.01, "npc"),
                cutoffLineType = "dotted", 
                cutoffLineCol = "#00000066",
                #titleLabSize = 10,
                #subtitleLabSize = 8,
                #captionLabSize = 8
)  + xlim(-.75, .75) 
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

##Barplot and Gene plot 5SE 5SD 5SB

table2 <- as.matrix(table(sce1@meta.data$condition, sce1@meta.data$cluster))
table2[1,] <- table2[1,] * table(sce1$condition)[2]/table(sce1$condition)[1]
table3 <- data.table(t(t(table2) / colSums(table2)))
colnames(table3) <- c("Condition", "Cluster", "Fraction")
ggplot(data=table3, aes(x=Cluster, y=Fraction, fill=Condition)) +
  geom_bar(stat="identity") + scale_fill_manual(values = c(ko_col, wt_col)) + ggtitle("Clusters by Condition")

sce1 <- sce2 
Idents(sce1) <- sce1$condition
sce1 <- subset(x = sce1, idents = "KO") 
th <- length(WhichCells(sce1, slot = 'counts', expression = TH > 0)) /  length(WhichCells(sce1, slot = 'counts'))
mki67 <- length(WhichCells(sce1, slot = 'counts', expression = MKI67 > 0)) /  length(WhichCells(sce1, slot = 'counts'))
map2 <- length(WhichCells(sce1, slot = 'counts', expression = MAP2 > 0)) /  length(WhichCells(sce1, slot = 'counts'))
all_ko <- c(map2, th, mki67)
sce1 <- sce2
Idents(sce1) <- sce1$condition
sce1 <- subset(x = sce1, idents = "WT")
th <- length(WhichCells(sce1, slot = 'counts', expression = TH > 0)) /  length(WhichCells(sce1, slot = 'counts'))
mki67 <- length(WhichCells(sce1, slot = 'counts', expression = MKI67 > 0)) /  length(WhichCells(sce1, slot = 'counts'))
map2 <- length(WhichCells(sce1, slot = 'counts', expression = MAP2 > 0)) /  length(WhichCells(sce1, slot = 'counts'))
all_wt <- c(map2, th, mki67)
temp <- data.table(
  Gene = c("MAP2", "TH", "MKI67"),
  Group =c("WT", "WT", "WT", "KO", "KO", "KO"),
  Values = c(all_wt, all_ko)
)

ggplot(temp, aes(x = Gene, y = Values, fill = Group)) + 
  geom_bar(stat = "identity", position = "dodge") + ggtitle("Fraction of Cells Expressing Each Gene") + ggeasy::easy_center_title() +
  theme(plot.title = element_text(face="bold")) +theme(text=element_text(size=8)) + scale_fill_manual(values = c(ko_col, wt_col))

##Supplement 5C

sce1 <- sce2
sce1.list = list(sce1[, sce1@meta.data$condition == "WT"], sce1[, sce1@meta.data$condition == "KO"])
f1 <-FeaturePlot(sce1.list[[1]], feature="MAP2", order=T, pt.size=1.0) + ggplot2::labs(title="MAP2 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f2 <-FeaturePlot(sce1.list[[2]], feature="MAP2", order=T, pt.size=1.0) + labs(title="MAP2 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#PBX1
f3 <-FeaturePlot(sce1.list[[1]], feature="PBX1", order=T, pt.size=1.0) + labs(title="PBX1 (WT)") + scale_colour_gradientn(colours =terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f4 <-FeaturePlot(sce1.list[[2]], feature="PBX1", order=T, pt.size=1.0) + labs(title="PBX1 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#MKI67
f5 <-FeaturePlot(sce1.list[[1]], feature="MKI67", order=T, pt.size=1.0) + labs(title="MKI67 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f6 <-FeaturePlot(sce1.list[[2]], feature="MKI67", order=T, pt.size=1.0) + labs(title="MKI67 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#SERT
f7 <-FeaturePlot(sce1.list[[1]], feature="SLC6A4", order=T, pt.size=1.0) + labs(title="SLC6A4 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f8 <-FeaturePlot(sce1.list[[2]], feature="SLC6A4", order=T, pt.size=1.0) + labs(title="SLC6A4 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#TPH1
f9 <-FeaturePlot(sce1.list[[1]], feature="TPH1", order=T, pt.size=1.0) + labs(title="TPH1 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f10 <-FeaturePlot(sce1.list[[2]], feature="TPH1", order=T, pt.size=1.0) + labs(title="TPH1 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
(f1/f2) | (f3/f4) | (f5/f6)

##FIGURE 5C

#load cellline data
cellline <- readRDS("cellline_clean.rds")
cellline <- as.Seurat(cellline)
## Warning: The column names of the counts matrix is NULL. Setting cell names to
## cell_columnidx (e.g 'cell_1').
## Warning: The column names of the data matrix is NULL. Setting cell names to
## cell_columnidx (e.g 'cell_1').
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from PC to PC_
## Warning: All keys should be one or more alphanumeric characters followed by an
## underscore '_', setting key to PC_
#La Manno dataset used a reference
combined_sce1 <- as.SingleCellExperiment(sce1)
combined_sce1$clusters <- sce1$clusters
combined_sce_lamanno <- scRNAseq::LaMannoBrainData("human-embryo")
## snapshotDate(): 2021-10-19
## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
## loading from cache
## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
## loading from cache
# SingleR() expects reference datasets to be normalized and log-transformed.
combined_sce_lamanno <- logNormCounts(combined_sce_lamanno) 
predictions2 <- SingleR(test=combined_sce1, assay.type.test=2, 
                        ref=combined_sce_lamanno, labels =combined_sce_lamanno$Cell_type, de.method="wilcox", clusters=combined_sce1$clusters)
predictions2$pruned.labels #Didn't assign
## [1] "hNbM"     "hNbM"     "hProgFPL" "hNbM"     "hProgFPL" "hProgFPL" "hPeric"
sce5 <- as.Seurat(combined_sce1)
Idents(sce5) <- sce5$clusters
sce5 <- RenameIdents(object = sce5, `1` = "hNbM", `2` = "hNbM", `3` = "hProgFPL", '4'= "hNbM", "5" = "hProgFPL", "6" = "hProgFPL", "7" = "hPeric")
DimPlot(sce5) + ggtitle("UMAP by Annotation")

combined_sce_lamanno1  <- subset(combined_sce_lamanno , ,Cell_type == c('hNbM') | Cell_type == c('hNbML1') | Cell_type == c('hProgFPL'))
predictions3 <- SingleR(test=combined_sce1, assay.type.test=2, 
                        ref=combined_sce_lamanno1, labels =combined_sce_lamanno1$Cell_type, de.method="wilcox")
combined_sce1$pruned.labels <- predictions3$pruned.labels
sce6 <- as.Seurat(combined_sce1)
v1 <- VlnPlot(sce6, "BAX", group.by = "pruned.labels", cols = c(wt_col, ko_col), pt.size = 0, split.by = "condition")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.
v2 <- VlnPlot(sce6, "TNFRSF12A", group.by = "pruned.labels",  cols = c(wt_col, ko_col),  pt.size = 0, split.by = "condition")
v3 <- VlnPlot(sce6, "JUN", group.by = "pruned.labels", cols = c(wt_col, ko_col),  pt.size = 0, split.by = "condition")
v1 + v2 + v3

cellline <- readRDS("cellline_clean.rds")
cellline <- as.Seurat(cellline)
## Warning: The column names of the counts matrix is NULL. Setting cell names to
## cell_columnidx (e.g 'cell_1').
## Warning: The column names of the data matrix is NULL. Setting cell names to
## cell_columnidx (e.g 'cell_1').
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from PC to PC_
## Warning: All keys should be one or more alphanumeric characters followed by an
## underscore '_', setting key to PC_
f1 <- FeaturePlot(sce1, feature="HES5", order=T, pt.size=1.0) + labs(title="HES5") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f1 <- FeaturePlot(sce1, feature="HES5", order=T, pt.size=1.0) + labs(title="HES5") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#ALTERNATE FIGURE REQUESTS
sce4 <- subset(sce1, ident == 1 |ident == 2 | ident == 4)
Idents(sce4) <-sce4$condition

markers2 <- FindAllMarkers(sce4, min.pct = 0.01, min.diff.pct = 0.01, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster KO
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
top2 <- markers2 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
  top2$avg_log2FC >= .10 & top2$p_val_adj < .01, '#00A08A',
  ifelse(top2$avg_log2FC <= -.10 & top2$p_val_adj < .01, '#F2AD00',
         'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
samplabels<- c("PEG10",
               "NMT1",
               "HSPA8", "RPS27L",
               "NDUFA12",
               "DYNLT1",
               "PDCD5",
               "PHPT1",
               "KRTCAP2", "SYT14",
               "SOX2",
               "HSPA1B",
               "HSPA1A","HSPA2A","DNAJB1","HSPA8")
EnhancedVolcano(top2,
                lab = top2$gene,
                x = 'avg_log2FC',
                y = 'p_val_adj',
                title = 'WT vs KO Top DE Genes',
                pCutoff = .01,
                FCcutoff = .10,
                pointSize = 1.0,
                labSize = 3.3,
                selectLab = samplabels,
                col=c('blue', 'green', 'black', 'red3', "yellow", "brown", "grey"),
                colAlpha = .8,
                colCustom = keyvals,
                drawConnectors = TRUE,
                widthConnectors = 0.1,
                lengthConnectors = unit(0.01, "npc"),
                cutoffLineType = "dotted", 
                cutoffLineCol = "#00000066",
                #titleLabSize = 10,
                #subtitleLabSize = 8,
                #captionLabSize = 8
)  + xlim(-.8, .8) 
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.